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A  Furiraff-TV1 subroutine  "is-grVen  which  may  be  used  to  generate 
random  observations  of  a  continuous  random  variable  from  a  table  of  dis¬ 
crete  observations.  This  subroutine  employs  Akima’s  algorithm  for  cubic 
spline  interpolation. 

The  compatibility  of  problem  and  algorithm  is  such  that  this  simple 
routine  gives  very  good  results. 


1,  The  Algorithm.  An  important  part  or  any 
simulation  study  is  the  sampling  of  a  ran¬ 
dom  variable  x  with  cummulative  distribu¬ 
tion  function  F (x) .  Often  F  is  specified 
only  by  a  table  of  (x,F(x)),  the  probability 
distribution  of  the  random  variable  y”F(x) 
is  the  uniform  distribution  U(0,1).  There 
are  satisfactory  pseudorandom  number  gen¬ 
erators  to  sample  from  U(0,1).  For  ex¬ 
ample  from  Knuth  [2] 


(1.1)  y 


■  314l5926yn+  27182818(mod  2J 


The  problem  of  samnling  from  F(x)  may  be 
looked  at,  then,  as  equivalent  to  finding 
an  inversion  formula 

(1.2)  x  *  F"*(y)  . 

To  approximate  F  *(y)  we  use  an  empiri¬ 
cally  obtained  table 

(1.3)  ((xl,F(xi)):x1  <x2  <  ...  <xfc) 

plus  a  good  interpolation  routine.  He 
should  have  in  (1.3)  that  F(x1)~0  and 
Ftx^)—!  .  In  this  way  it  13  possible  to 
Lake  random  values  from  U(0,1)  and  obtain 
their  approximate  Inverses  as  random  samples 
from  "nearly  F(x)." 


STJSSaTia) »  lml^1  ~~ 


The  sampling  problem  now  reduces  to  one 
of  interpolation  subject  to  a  montonicity 
constraint  or  "isotone  interpolation." 

Neither  the  standard  polynomial  interpolation 
nor  the  cubic  spline  interpolation  have  a 
chance  of  accomplishing  this  objective.  For 
this  reason  (and  others)  the  standard  approach 
in  the  literature  is  piecewise  linear  inter¬ 
polation.  Observe  that  piecewise  liner  inter¬ 
polation  maintains  monotonicity  but,  in 
general,  at  the  expense  of  accuracy. 

It  is  the  purpose  of  this  note  merely  to 
point  out  that  Akiaa's  {1J  quasi-Hermite 
piecewise  cubic  interpolation  should  be  very 
effective  on  this  problem;  since  it  was  de¬ 
signed  to  handle  similar  problems.  A  random 
number  generator  (RNG)  is  advocated  which 
uses  Akima  interpolation.  The  user  reads  in 
a  table 

((*(1),F(*(i))):x(0)  <x(1)  <  ...  <x(k)  <x(fcfl)). 

The  routine  then  computes  and  stores  the  co¬ 
efficients  of  x  =  F*f(y).  Calls  on  the  ran¬ 
dom  number  generator  will  first  obtain  a 
pseudorandom  number  v  from  U(0,1).  If 
y<F(x(1)),  the  random  number  given  is  , 

if  y>F(x^),  the  random  number  given  is  x^  . 
If  P(*(i)  <y<F(x.  ),  the  random  number  given 
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isF  6r) .  A  listing  of  the  random  number  gen¬ 
erator  is  given  in  Appendix  II. 

Naturally,  the  quality  of  the  random 
number  generator  proposed  is  determined  by 
the  ability  of  F"1  to  approximate  a  var¬ 
iety  of  cumulative  distribution  functions. 

In  Appendix  I,  we  show  how  the  interpola¬ 
tion  errors  of  F*^(y)  compare  with  those 
of  linear  interpolation  in  the  case  of  the 
standard  normal,  Cauchy  and  double  exponen¬ 
tial  distributions  with  location  0  and 
scale  1.  Exact  values  of  F-l(y)  are 
assumed  given  for  y  between  .50  and  .98 
in  increments  of  .01  and  for  y  between 
.98  and  .995  in  increments  of  .005.  In¬ 
terpolated  values  of  F'*  and  linear  inter¬ 
polation  are  obtained  at  points  between  the 
mesh  values  and  the  relative  errors  com¬ 
puted  and  compared. 
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APPENDIX  I 
NORMAL  DISTRIBUTION 

F  (x)*=  x»  E(AKIMAS)=  E  (LINEAR) 


0.50333 

0.00833 

0.52333 

0.05837 

0.54333 

0.10858 

0.56333 

0.15910 

0.58333 

0.21005 

0.60333 

0.26157 

0.62333 

0.31381 

C. 64333 

0.36694 

0.66333 

0.42114 

0.68333 

0.47662 

0.70333 

0.53362 

0.72333 

0.59241 

0.74333 

0.65334 

0.76333 

0.71680 

0.78333 

0.78329 

0.80333 

0.85343 

0.82333 

0.92805 

0.84333 

1.00823 

0.86333 

1.09546 

0.88333 

1.19193 

0.90333 

1.30097 

0.?2 .ii 

• .47812 

0.94333 

1.58372 

0.96333 

1.79115 

0.98333 

2.12349 

0.99333 

2.47156 

0.00006 

0.00011 

0.00003 

0.00020 

0.00002 

0.00033 

0.00006 

0.00044 

0.00004 

0.00061 

0.00001 

0.00080 

0.00001 

0.00096 

0.00000 

0.00114 

0.00004 

0.00129 

0.00005 

0.00149 

0.00004 

0.00172 

0.00002 

0.00200 

0.00000 

0.00231 

0.00007 

0.00258 

0.00004 

0.00300 

0.00005 

0.00346 

0.00007 

0.00401 

0.00012 

0.00469 

0.00011 

0.00561 

0.00018 

0.00679 

0.00021 

0.00853 

0.00030 

0.01115 

0.00034 

0.01569 

0.00185 

0.02539 

0.00336 

0. 02 805 

0.02432 

0.07186 

F(x)  = 

x= 

E(AKIMAS)= 

E  (LINEAR): 

0.50333 

0.01047 

0.00002 

0.00010 

0.52333 

0.07344 

0.00005 

0.00054 

0.54333 

0.13698 

0.00005 

0.00098 

0.56333 

0.20164 

0.00005 

0.00143 

0.58333 

0.26795 

0.00005 

0.00190 

0.60333 

0.33654 

0.00006 

0.00237 

0.62333 

0.40809 

0.00006 

0.00287 

0.64333 

0.48342 

0.00007 

0.00340 

0.66333 

0.56347 

0.00008 

0.00396 

0.68333 

0.64941 

0.00009 

0.00456 

0.70333 

0.74267 

0.00011 

0.00521 

0.72333 

0.84507 

0.00013 

0.00592 

0.74333 

0.95897 

0.00015 

0.00672 

0.76333 

1.08749 

0.00017 

0.00761 

0.78333 

1.23490 

0.00021 

0.00865 

0.80333 

1.40714 

0.00025 

0.00985 

0.82333 

1.61283 

0.00030 

0.01128 

0.84333 

1.86499 

0.00037 

0.01304 

0.86333 

2.18419 

0.00046 

0.01526 

0.88333 

2.60509 

0.00057 

0.01820 

0.90333 

3.19100 

0.00068 

0.02229 

0.92333 

4.07127 

0.00072 

0.02844 

0.94333 

5.55776 

0.00021 

0.03882 

0.96333 

8.64274 

0.00113 

0.06036 

0.98333 

19.08113 

0.02053 

0.06658 

0.99333 

47.73947 

0.08864 

0.16661 

DOUBLE  EXP 

DISTRIBUTION 

F(x)  = 

X“ 

E(AKIMAS)- 

E(  LINEAR) 

0.50333 

2.09951 

0.00002 

0.00223 

0.52333 

2.22281 

0.00003 

0.00233 

0.54333 

2.35140 

0.00002 

0.00244 

0.56333 

2.48575 

0.00002 

0.00255 

0.58333 

2.62641 

0.00002 

0.00268 

0.60333 

2.77398 

0.00003 

0.00281 

0.62333 

2.92918 

0.00003 

0.00295 

0.64333 

3.09286 

0.00004 

0.00312 

0.66333 

3.26599 

0.00003 

0.00331 

0.68333 

3.44972 

0.00004 

0.00351 

0.70333 

3.6'544 

0.00005 

0.00375 

0.72333 

3.85483 

0.00006 

0.00402 

0.74333 

4.07993 

0.00006 

0.00434 

0.76333 

4.32331 

0.00007 

0.00471 

0.78333 

4.58819 

0.00009 

0.00514 

0.80333 

4.87874 

0.00010 

0.00567 

0.82333 

5.20047 

0.00013 

0.00631 

0.84333 

5.56091 

0.00015 

0.00712 

0.86333 

5.97063 

0.00019 

0.00817 

0.88333 

6.44530 

0.00024 

0.00957 

0.90333 

7.00946 

0.00031 

0.011 

57 

0.92333 

7.70487 

0.00040 

0.01461 

0.94333 

8.61171 

0.00042 

0.01983 

0.96333 

9.91767 

0.00215 

0.03088 

0.98333 

12.28305 

0.00450 

0.03291 

0.99333 

is  nrtos 

n  miM 

- - 
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APPENDIX  II 


_C  INTERPOL  AT  J  pN  O.F_THE  INVERSE  NORML  O I  STM  I  OUT  I  ON 
COMMON  N  .P~t  (SO)  ,P2(  SO)  *P3<  50)  .P4|50>  .XI  50  ) 
c  N  IS  The  NUMBER  OK  INTERPOLATING  KNOTS,  n  should  oe  less  or  CO  TO 
N  m  50 
CALL  I N|T 

C  M  IS  THE  NUMOER  OF  RANDOM  NUMBERS  TO  OE  GENERATED  M  SMALLER  THAN 
M  *  20 

'call  rngcmi 

CALL  EXIT 
ENO 

SU8.R0UT  |  NE  IN  IT 

C  •** *•••••»» *  rY*  *•  *" 

C  THIS  SUBROUTINE  GENERATES  THE  PARAMETERS  OR  T«£  CUB ICS  THAT 
C  INTERPOLATE  THE  DATA  KNOTS 

C  ARRAYS  PI.PP.P3.P4  CONTAIN  THF  PARAMETERS 

COMMON  N.P1  <50)  .!>?(  SO)  .P3ISOI  .PA(SO)  «A(S0| 

100  FORMAT  (2F12.6) 

200  FORMAT  (2X.10FI2.S) 

C  •**•**•  READ  DATA  ****•»*»•••*•*••••*«••*•• 

READ  too. ( (x« t) ,P4( I) I .1*1 ,N) 

C  ••***•*»**  READ  OATA  FINISHES  ••**•*•••«• 

N2  =  N-? 

3 1*1  PA  (  2  )  -P4  M  l  )  ✓  (  x  <  2  )-X<  1)> 

S2«<  PA(  3  )-PA  (  2)  |  /(  Y<  3)  -X  (Yf) 

S3=CPA( 3 )-PA( 41 )/(X(3)-X(4) ) 

DO  10  1=3. N2 

S4  *  (P4(I.?>  -P4C  I.t))/(X(  1421  -  X(IAU) 

Ml  *  A0SIS2-S1) 

_W3  *_  AS  S  (  54 -S  3  > 

IF  ~(MI  *W3  >  20  *  3 0  «  20 
20  P3III  =  (W3*S2  ♦  MI AS3I/(M|+W3) 

CO  TO  40 

30  P3II)  =  • 5*1 S2+S3 ) 

40  SI  *  S2 

S2  _■  S3 
10  S3  *  S4 

00  60  I=3.N2 
A*-  X  t  I  ♦  I  1  -  X<  I  ) 

p|<if  *  P3<i*  n  -  2.0*1  pa  f  r  *  i  /  -  **4(  r  jisai/caaai 

P2<  I  la(3.0*(P4(  I  +  l  *-P4(  I  )  )  /A  -  2  •  0*P3f  I  )  -  P3(I*11)/A 

_$* 3.  continue 

RETURN 
END 

SUBROUTINE  RNG(M). 

COMMON  N.Pt (50) ,02(50) .P3< 50 ) .PA C 50 ) *X(50) 

REALAS  RN(?0).I3A.I35 
100  FORMAT  (16. OF  12.6) 

C  AAA  AAA  RANDOM  NUMBER  GENERATOR  AAAA**AA 
13A  *  57721566 
.13®  J«_2.0AA35 
od*’io“i*i  *M 

I  3 A  *  (  3 1 A  1 5926  • 1 3A  ♦  27182628  ) 

134  «  (DMOOC I  34. 135) ) 

RN( I )  *  134/135 
PRINT  lOO.I.RN(I) 

IF  CRNltJ .LE.XC3))  RN( I >  •  X(3I 
«0  IF  (RNlVl  .GE«X<N-2) )  RN(  II*  X«N-2) 

•  IC  *  1 
NS  *  N-3 
OO  20  I  -3.N3 
DO  30  4*1  .  M 
.V  »  RNIJ)  -Mill 
*  IF  (RN(J)-X(I))  40.50.50 
SO  IF  <RN(J)  -  XdAlll  60.60.40 

60  VRN  ■  P I (  I ) A V* y* V  ♦  P2 (  I  I  A V AV  ♦  P3(I|AV  ♦  P4f I  1 
C  VRN  IS  the  GENERATED  RANDOM  OBSERVATION 
PRINT  100.IC.RNC Jl.TRN 
IC  •  IC  ♦  1 


40 

CONTINUE 

tP  UC.GT 

30 

COHf INUE 

to 

CONT | NUC 

0 

CONTINUE 

RETURN 

ENO 


/ 
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